function S = slpca(X, varargin)
%SLPCA Learns a PCA model from training samples
%
% [ Syntax ]
%   - S = slpca(X)
%   - S = slpca(X, ...)
%
% [ Arguments ]
%   - X:        the training sample matrix
%   - S:        the struct representing the learned PCA
%               fields:
%               \{:
%                  - sampledim:     the dimension of samples
%                  - feadim:        the dimension of subspace
%                                   feature
%                  - vmean:         the mean vector
%                  - P:             the projection matrix
%                                   (preserved eigenvectors)
%                  - eigvals:       the vector of preseved eigen values
%                  - total:         the total energy
%                  - preserved:     the preserved energy
%                  - residue:       the residue energy
%                  - energyratio:   the ratio of preserved energy
%               \:}
%
% [ Description ]
%   - S = slpca(X) learns a PCA model from the samples X with default options.
%
%   - S = slpca(X, ...) learns a PCA model from the samples X according to
%     the properties specified:
%     \{:
%         - method:     The method used in training the PCA model. 
%                       \{:
%                           - auto:  automatically select the best method
%                           - std:   use standard way based on covariance
%                           - svd:   use SVD to decompose sample matrix
%                           - trans: use a transposed way, which solves the
%                                    eigen-problem of X^T * X, and then 
%                                    transform back to the original domain                           
%                       \:}
%                       By default, method is 'auto'.
%
%         - preserve:   The cell array in form of {cri} or {cri, v} to specify
%                       which eigenvalues (and corresponding eigenvectors) are
%                       preserved.
%                       Here, cri is the name of the criteria as listed below, 
%                       v is the parameter value.
%                       \{:
%                           - rank:   preserve all eigenvalues not smaller 
%                                     than eps(max(eigvals))
%                           - num:    preserve the v leading eigenvalues
%                           - ratio:  preserve all eigenvalues not smaller
%                                     than (v * max(eigvals))
%                           - energy: preserve energy of ratio not less
%                                     than v.                           
%                       \:}
%                       By default, preserve is {'rank'}.
%
%         - selfun:     The function handle to select preserved eigenvalues.
%                       It should support the syntax as
%                          $ inds = f(eigvals) $
%                       It outputs the indices of the selected eigenvalues
%                       when input all the eigenvalues.
%                       default = [].
%
%         - weights:    The weights of the samples. default = [].
%
%         - vmean:      The pre-computed mean vector of the samples. 
%                       default = []. It can be specified as a column
%                       vector or zero which means that all samples are
%                       centralized in advance.                       
%     \:}
%
% [ Remarks ]
%   - If selfun is specified as a function handle, then the function will
%     use selfun to select the preserved eigenvalues, otherwise the
%     function uses preserve option to select eigenvalues.
%                           
% [ History ]
%   - Created by Dahua Lin on Apr 24, 2006
%   - Modified by Dahua Lin on Apr 25, 2006
%       - Extract the dimension determination part to an independent
%         function, in order to offer more flexible preservation process,
%         and make it reusable in more functions.
%   - Modified by Dahua Lin on Aug 17, 2006
%       - Add a field energyratio
%   - Modified by Dahua Lin, on Sep 10, 2006
%       - replace sladd by sladdvec to increase efficiency
%   - Modified by Dahua Lin, on Oct 1, 2006
%       - add the support of input vmean in various forms
%   - Modified by Dahua Lin, on Jul 16, 2007
%       - refactor the code and document
%       - add support to user-supplied eigenvalue selection function
%       - add more fields related to energy preservation to the output
%         struct
%

%% parse and verify input arguments

% X
assert(isfloat(X) && ndims(X) == 2, 'sltoolbox:slpca:invalidarg', ...
    'X should be a 2D float/double matrix.');
[d, n] = size(X);


% for options
opts = struct( ...
    'method', 'auto', ...
    'preserve', {{'rank'}}, ...
    'selfun', [], ...
    'weights', [], ...
    'vmean', []);
if ~isempty(varargin)
    opts = setopts(opts, varargin{:});                
end

% check options

method_names = {'auto', 'std', 'svd', 'trans'};
method_funcs = {@pca_auto, @pca_std, @pca_svd, @pca_trans};

assert(ischar(opts.method), 'sltoolbox:slpca:invalidopt', ...
    'The option method should be a string.');
idx_method = find(strcmp(opts.method, method_names));
assert(~isempty(idx_method), 'The option method can not be set to %s', opts.method);

fh_learnpca = method_funcs{idx_method};

if isempty(opts.selfun)
    assert(iscell(opts.preserve) && (numel(opts.preserve) == 1 || numel(opts.preserve) == 2), ...
        'sltoolbox:slpca:invalidopt', 'The option preserve should be a cell array of length 1 or 2.');
    fh_seleig = geteigpreserve(opts.preserve{:});
    
else
    assert(isa(opts.selfun, 'function_handle'), 'sltoolbox:slpca:invalidopt', ...
        'The option selfun should be either empty or a function handle.');
    
    fh_seleig = opts.selfun;
end

if ~isempty(opts.weights)
    assert(isnumeric(opts.weights) && isvector(opts.weights), ...
        'sltoolbox:slpca:invalidopt', 'The option weights should be either empty or a numeric vector.');
    assert(length(opts.weights) == n, 'sltoolbox:slpca:invalidopt', ...
        'sltoolbox:slpca:invalidopt', 'The length of weights should equal the number of samples.');
    assert(all(opts.weights >= 0) && any(opts.weights > 0), ...
        'sltoolbox:slpca:invalidopt', 'All weights should be nonnegative and some weights positive.');
end

if ~isempty(opts.vmean) && ~isequal(opts.vmean, 0)
    assert(isnumeric(opts.vmean) && isequal(size(opts.vmean), [d 1]), ...
        'sltoolbox:slpca:invalidopt', ...
        'The option vmean should be either empty, zero or a d x 1 column vector');        
end


%% main

% data centralization
if isempty(opts.vmean)
    vmean = slmean(X, opts.weights);
else
    vmean = opts.vmean;
end

if ~all(vmean == 0)
    X = bsxfun(@minus, X, vmean);
end

% learn full size PCA

if isempty(opts.weights)
    w = opts.weights;
else
    w = reshape(opts.weights, [1 n]);
    w = w / sum(w);
end

[P, evals] = fh_learnpca(X, w);

kmax = min(d, n-1);
if kmax < size(P, 2)
    P = P(:, 1:kmax);
    evals = evals(1:kmax);
end
evals(evals < 0) = 0;
total_energy = sum(evals);

% preserve principal components

inds = fh_seleig(evals);
evals = evals(inds);
P = P(:, inds);

S = makepcamodel(vmean, evals, P, total_energy);



%% Sub functions for Learning PCA from centralized samples and normalized weights

function [P, evals] = pca_auto(X, w)

[d, n] = size(X);
if d <= n
    [P, evals] = pca_std(X, w);
else
    [P, evals] = pca_trans(X, w);
end


function [P, evals] = pca_std(X, w)

if isempty(w)
    C = X * X' / size(X, 2);
else
    C = bsxfun(@times, X, w) * X';
end 
[evals, P] = slsymeig(C);


function [P, evals] = pca_svd(X, w)

if isempty(w)
    [P, D] = svd(X, 0);
    evals = diag(D) .^ 2 / size(X, 2);
else
    [P, D] = svd(bsxfun(@times, X, sqrt(w)), 0);
    evals = diag(D) .^ 2;
end    


function [P, evals] = pca_trans(X, w)

if isempty(w)
    Ct = X' * X / size(X, 2);
else
    Z = bsxfun(@times, X, sqrt(w));
    Ct = Z' * Z;
end
[evals, P] = slsymeig(Ct);

if isempty(w)
    P = slnrmvecs(X * P);
else
    P = slnrmvecs(Z * P);
end
    
    




